home *** CD-ROM | disk | FTP | other *** search
/ Amiga Format CD 46 / Amiga Format CD46 (1999-10-20)(Future Publishing)(GB)[!][issue 1999-12].iso / -in_the_mag- / reader_requests / scilab / tests / matode.dia.ref < prev    next >
Text File  |  1999-09-16  |  11KB  |  590 lines
  1.  
  2. Leps=1.e-6;
  3.  
  4. //     dy1/dt = -0.04*y1 + 1.e4*y2*y3
  5.  
  6. //     dy2/dt =0.04*y1 - 1.e4*y2*y3 - 3.e7*y2**2
  7.  
  8. //     dy3/dt = 3.e7*y2**2
  9.  
  10. // on the interval from t = 0.0 to t = 4.e10, with initial conditions
  11.  
  12. // y1 = 1.0, y2 = y3 = 0.  the problem is stiff.
  13.  
  14. //     definition of rhs
  15.  
  16. deff('[ydot]=f(t,y)',...
  17.    'f1=-0.04*y(1)+1e4*y(2)*y(3),...
  18.     f3=3e7*y(2)*y(2),...
  19.     ydot=[f1;-f1-f3;f3]')
  20.  
  21. //     jacobian
  22.  
  23. deff('[jac]=j(t,y)',...
  24. 'jac(1,1)=-0.04;jac(1,2)=1.e4*y(3);jac(1,3)=1.e4*y(2),...
  25.  jac(3,1)=0;jac(3,2)=6.e7*y(2);jac(3,3)=0;...
  26.  jac(2,1)=-jac(1,1);jac(2,2)=-jac(1,2)-jac(3,2);jac(2,3)=-jac(1,3);')
  27.  
  28. //
  29.  
  30. y0=[1;0;0];t0=0;t1=[0.4,4];nt=size(t1,'*');
  31.  
  32. //    solution
  33.  
  34. yref=[0.9851721 0.9055180;0.0000339 0.0000224;0.0147940 0.0944596];
  35.  
  36. //
  37.  
  38. //  1. fortran called by fydot, without jacobian
  39.  
  40. y1=ode(y0,t0,t1,'fex');
  41.  
  42. if maxi(y1-yref) > Leps then bugmes();quit;end
  43.  
  44. //  2. fortran called by fydot, type given (stiff), no jacobian
  45.  
  46. y2=ode('stiff',y0,t0,t1,'fex');
  47.  
  48. if maxi(y2-yref) > Leps then bugmes();quit;end
  49.  
  50. //  3. fortran called by fydot , fjac, type given
  51.  
  52. y3=ode('stiff',y0,t0,t1,'fex','jex');
  53.  
  54. if maxi(y3-yref) > Leps then bugmes();quit;end
  55.  
  56. //   hot restart
  57.  
  58. [z,w,iw]=ode('stiff',y0,0,0.4,'fex','jex');
  59.  
  60. z=ode('stiff',z,0.4,4,'fex','jex',w,iw);
  61.  
  62. if maxi(z-y3(:,2)) > %eps then bugmes();quit;end
  63.  
  64.  
  65. [y1,w,iw]=ode(y0,t0,t1(1),'fex');
  66.  
  67. y2=ode(y0,t1(1),t1(2:nt),'fex',w,iw);
  68.  
  69. if maxi([y1 y2]-yref) > Leps then bugmes();quit;end
  70.  
  71.  
  72. [y1,w,iw]=ode(y0,t0,t1(1),'fex','jex');
  73.  
  74. y2=ode(y0,t1(1),t1(2:nt),'fex','jex',w,iw);
  75.  
  76. if maxi([y1 y2]-yref) > Leps then bugmes();quit;end
  77.  
  78.  
  79. //   variation of tolerances
  80.  
  81. atol=[0.001,0.0001,0.001];rtol=atol;
  82.  
  83. //    externals
  84.  
  85. comp(f),comp(j)
  86.  
  87. //  4. type given , scilab lhs ,jacobian not passed
  88.  
  89. y4=ode('stiff',y0,t0,t1(1),atol,rtol,f);
  90.  
  91. if maxi(y4(:,1)-yref(:,1)) > 0.01 then bugmes();quit;end
  92.  
  93. //  5. type non given, rhs and scilab jacobian
  94.  
  95. y5=ode(y0,t0,t1,f,j);
  96.  
  97. if maxi(y5-yref) > Leps then bugmes();quit;end
  98.  
  99. //  6. type given (stiff),rhs and jacobian  by scilab
  100.  
  101. y6=ode('stiff',y0,t0,t1,0.00001,0.00001,f,j);
  102.  
  103. if (y6-yref) > 2*0.00001 then bugmes();quit;end
  104.  
  105. //  7. matrix rhs, type given(adams),jacobian not passed
  106.  
  107. //
  108.  
  109. a=rand(3,3);ea=exp(a);
  110.  
  111. deff('[ydot]=f(t,y)','ydot=a*y')
  112. Warning :redefining function: f                       
  113.  
  114.  
  115. comp(f);
  116.  
  117. t1=1;y=ode('adams',eye(a),t0,t1,f);
  118.  
  119. if maxi(ea-y) > Leps then bugmes();quit;end
  120.  
  121. //
  122.  
  123. //   DAE's
  124.  
  125. //     dy1/dt = -0.04*y1 + 1.e4*y2*y3
  126.  
  127. //     dy2/dt =0.04*y1 - 1.e4*y2*y3 - 3.e7*y2**2
  128.  
  129. //       0.   = y1 + y2 + y3 - 1.
  130.  
  131. //  y1(0) = 1.0, y2(0) = y3(0) = 0.
  132.  
  133. // scilab macros
  134.  
  135. deff('[r]=resid(t,y,s)','...
  136. r(1)=-0.04*y(1)+1.d4*y(2)*y(3)-s(1),...
  137. r(2)=0.04*y(1)-1.d4*y(2)*y(3)-3.d7*y(2)*y(2)-s(2),...
  138. r(3)=y(1)+y(2)+y(3)-1')
  139.  
  140. deff('[p]=aplusp(t,y,p)','...
  141. p(1,1)=p(1,1)+1,...
  142. p(2,2)=p(2,2)+1')
  143.  
  144. deff('[p]=dgbydy(t,y,s)','[-0.04,1.d4*y(3),1.d4*y(2);...
  145. 0.04,-1.d4*y(3)-6.d7*y(2),-1.d4*y(2);...
  146. 1,1,1]')
  147.  
  148. comp(resid);comp(aplusp);comp(dgbydy);
  149.  
  150. //        %ODEOPTIONS tests
  151.  
  152. //
  153.  
  154. rand('seed',0);rand('normal');
  155.  
  156. nx=20;A=rand(nx,nx);A=A-4.5*eye;
  157.  
  158. deff('y=f(t,x)','y=A*x')
  159. Warning :redefining function: f                       
  160.  
  161.  
  162. deff('J=j(t,x)','J=A')
  163. Warning :redefining function: j                       
  164.  
  165.  
  166. x0=ones(nx,1);t0=0;t=[1,2,3,4,5];
  167.  
  168. nt=size(t,'*');
  169.  
  170. eAt=exp(A*t(nt));
  171.  
  172.  
  173. //        Test itask=%ODEOPTIONS(1)
  174.  
  175.  
  176. //itask=1  --->  usual call (t=vector of instants, solution at all t)
  177.  
  178. //========================
  179.  
  180. itask=1;tcrit=0;h0=0;hmax=0;hmin=0;ixpr=0;mxstep=0;maxordn=0;maxords=0;
  181.  
  182. jacflag=2;ml=-1;mu=-1;
  183.  
  184. %ODEOPTIONS=[itask,tcrit,h0,hmax,hmin,jacflag,mxstep,maxordn,maxords,ixpr,ml,mu];
  185.  
  186. xf=ode(x0,t0,t,f);   //lsoda
  187.  
  188. if norm(xf(:,nt)-eAt*x0)>Leps then bugmes();quit;end
  189.  
  190. xfj=ode(x0,t0,t,f,j);   //lsoda with jacobian
  191.  Warning: Jacobian external is given, but 
  192.  not used!,  see %ODEOPTIONS(6)
  193.  
  194.  
  195. if norm(xfj(:,nt)-eAt*x0)>Leps then bugmes();quit;end
  196.  
  197.  
  198.  
  199. //itask=2;   --->  solution at mesh points  ---> t=tmax
  200.  
  201. //========================
  202.  
  203. %ODEOPTIONS(1)=2;tmax=t(5);
  204.  
  205. xft=ode(x0,t0,tmax,f);
  206.  
  207. [p,q]=size(xft);
  208.  
  209. xlast=xft(2:nx+1,q);
  210.  
  211. if xft(1,q)<tmax then bugmes();quit;end
  212.  
  213. if norm(xlast-exp(A*xft(1,q))*x0)>Leps then bugmes();quit;end
  214.  
  215.  
  216. //itask=3;   ---> solution at first mesh point beyond t=tmax
  217.  
  218. %ODEOPTIONS(1)=3;
  219.  
  220. x3=ode(x0,t0,tmax,f);
  221.  
  222. if norm(x3(2:nx+1)-xlast,1)>Leps then bugmes();quit;end
  223.  
  224.  
  225. //itask=4; test with %ODEOPTIONS(2)=tcrit
  226.  
  227. %ODEOPTIONS(1)=4; //---> computation at all t and t>tcrit are not called
  228.  
  229. tcrit=2.5;%ODEOPTIONS(2)=tcrit;
  230.  
  231. chk=0;
  232.  
  233. deff('y=fcrit(t,x)',['if t<=tcrit then'
  234.                       ' y=A*x;'
  235.                       'else'
  236.                       ' y=A*x;chk=resume(1);end'])
  237.  
  238. x42=ode(x0,t0,t,fcrit);
  239.  Warning: integration up to tcrit
  240.  
  241.  
  242. if chk==1 then bugmes();quit;end
  243.  
  244. [p,q]=size(x42);
  245.  
  246. if norm(x42(:,q)-ode(x0,t0,tcrit,f),1)>Leps then bugmes();quit;end
  247.  
  248.  
  249. //itask=5; test with %ODEOPTIONS(2)=tcrit
  250.  
  251. %ODEOPTIONS(1)=5;  //---> computation at mesh points and t>tcrit are not called
  252.  
  253. %ODEOPTIONS(6)=2;  // Estimated jacobian
  254.  
  255. chk=0;
  256.  
  257. x52=ode(x0,t0,2.3,fcrit);
  258.  
  259. if chk==1 then bugmes();quit;end
  260.  
  261. [p,q]=size(x52);
  262.  
  263. if x52(1,q)>tcrit then bugmes();quit;end
  264.  
  265.  
  266. //test of %ODEOPTIONS(3:5)=[h0,hmax,hmin]
  267.  
  268. %ODEOPTIONS(1)=1;
  269.  
  270. h0=0.0;hmax=0.1;hmin=0.0001;
  271.  
  272. %ODEOPTIONS(3:5)=[h0,hmax,hmin];
  273.  
  274. x35=ode(x0,t0,t,f);
  275.  
  276. if norm(x35-xf,1)>Leps then bugmes();quit;end
  277.  
  278.  
  279. //test of %ODEOPTIONS(6)=jacflag
  280.  
  281. %ODEOPTIONS(6)=1;//Jacobian given
  282.  
  283. %ODEOPTIONS(3:5)=[0 0 0];
  284.  
  285. x61=ode('st',x0,t0,t,f,j);   //with Jacobian
  286.  
  287. if norm(x61-xf,1)>10*Leps then bugmes();quit;end
  288.  
  289.  
  290.  
  291. %ODEOPTIONS(6)=0; // jacobian nor called nor estimated
  292.  
  293. x60=ode('st',x0,t0,t,f,j);   //Jacobian not used (warning)
  294.  Warning: Jacobian external is given, but 
  295.  not used!,  see %ODEOPTIONS(6)
  296.  
  297.  
  298. x60=ode('st',x0,t0,t,f);    //Jacobian not used
  299.  
  300. if norm(x60-x61,1)>10*Leps then bugmes();quit;end
  301.  
  302.  
  303.  
  304. %ODEOPTIONS(6)=1;//Jacobian estimated
  305.  
  306. x60=ode('st',x0,t0,t,f)  ;
  307.  Warning: No Jacobian external given but 
  308.  one is required by %ODEOPTIONS(6) value!
  309.  
  310.  
  311. if norm(x60-x61,1)>10*Leps then bugmes();quit;end
  312.  
  313.  
  314. //test of %ODEOPTIONS(6)=jacflag   (adams)
  315.  
  316. %ODEOPTIONS(6)=1;//with given Jacobian
  317.  
  318. x60=ode('ad',x0,t0,t,f,j) ;
  319.  
  320. if norm(x60-x61,1)>10*Leps then bugmes();quit;end
  321.  
  322.  
  323.  
  324. %ODEOPTIONS(6)=0;// jacobian nor called nor estimated
  325.  
  326. x60=ode('ad',x0,t0,t,f,j);   //Jacobian not used (warning)
  327.  Warning: Jacobian external is given, but 
  328.  not used!,  see %ODEOPTIONS(6)
  329.  
  330.  
  331. x60=ode('ad',x0,t0,t,f);    //Jacobian not used
  332.  
  333. if norm(x60-x61,1)>10*Leps then bugmes();quit;end
  334.  
  335.  
  336. // test lsoda
  337.  
  338. %ODEOPTIONS(6)=2;// jacobian  estimated
  339.  
  340. x60=ode(x0,t0,t,f);
  341.  
  342. if norm(x60-x61,1)>10*Leps then bugmes();quit;end
  343.  
  344.  
  345. %ODEOPTIONS(6)=1;//Jacobian estimated
  346.  
  347. x60=ode(x0,t0,t,f);
  348.  Warning: No Jacobian external given but 
  349.  one is required by %ODEOPTIONS(6) value!
  350.  
  351.  
  352. if norm(x60-x61,1)>10*Leps then bugmes();quit;end
  353.  
  354.  
  355.  
  356. //   Banded Jacobian
  357.  
  358. itask=1;tcrit=0;h0=0;hmax=0;hmin=0;ixpr=0;mxstep=0;maxordn=0;maxords=0;
  359.  
  360. //provisional values as default
  361.  
  362. jacflag=2;ml=-1;mu=-1;
  363.  
  364.  
  365. %ODEOPTIONS=[itask,tcrit,h0,hmax,hmin,jacflag,mxstep,maxordn,maxords,ixpr,ml,mu];
  366.  
  367. jacflag=2;%ODEOPTIONS(6)=jacflag;   //Banded Jacobian, given
  368.  
  369. nx=20;A=diag(-[1:nx]);x0=ones(nx,1);
  370.  
  371. t0=0;t=[1,2,3,4,5];
  372.  
  373. for k=1:nx-1, A(k,k+1)=1;end
  374.  
  375. for k=1:nx-2, A(k,k+2)=-2;end
  376.  
  377. for k=1:nx-1, A(k+1,k)=-1;end
  378.  
  379. for k=1:nx-2, A(k+2,k)=2;end
  380.  
  381. for k=1:nx-3, A(k+3,k)=-3;end
  382.  
  383. deff('xd=f(t,x)','xd=A*x')
  384. Warning :redefining function: f                       
  385.  
  386.  
  387. ml=3;mu=2;
  388.  
  389. %ODEOPTIONS(11:12)=[ml,mu];
  390.  
  391. for i=1:nx;
  392.     for j=1:nx;
  393. if A(i,j)<>0 then J(i-j+mu+1,j)=A(i,j);end
  394. end;end;
  395. Warning :redefining function: j                       
  396.  
  397.  
  398. // J is a ml+mu+1 x ny matrix.
  399.  
  400. // Column 1 of J is made of mu zeros followed by df1/dx1, df2/dx1, df3/dx1,
  401.  
  402. // i.e. 1 + ml possibly nonzeros entries.
  403.  
  404. // Column 2 of J is made of mu-1 zeros followed by df1/dx2, df2/dx2,0...
  405.  
  406. // etc...
  407.  
  408. %ODEOPTIONS(6)=1;%ODEOPTIONS(11:12)=[-1,-1];
  409.  
  410. deff('jj=j1(t,x)','jj=A')
  411.  
  412. xnotband=ode('st',x0,t0,t,f,j1);
  413.  
  414.  
  415.  
  416. %ODEOPTIONS(6)=4;//banded jacobian external given
  417.  
  418. %ODEOPTIONS(11:12)=[3,2];
  419.  
  420. deff('jj=j(t,x)','jj=J')
  421.  
  422. xband=ode('st',x0,t0,t,f,j);
  423.  
  424. if norm(xnotband-xband,1)>Leps then bugmes();quit;end
  425.  
  426.  
  427. %ODEOPTIONS(6)=5;//banded jacobian evaluated
  428.  
  429. %ODEOPTIONS(11:12)=[3,2];
  430.  
  431. deff('jj=j(t,x)','jj=J')
  432.  
  433. xband=ode('st',x0,t0,t,f,j);
  434.  Warning: Jacobian external is given, but 
  435.  not used!,  see %ODEOPTIONS(6)
  436.  
  437.  
  438. if norm(xnotband-xband,1)>Leps then bugmes();quit;end
  439.  
  440.  
  441.  
  442. //            Test of %ODEOPTIONS(7)
  443.  
  444. //%ODEOPTIONS(7)=mxstep  ---> maximum number od steps allowed
  445.  
  446. itask=1;tcrit=0;h0=0;hmax=0;hmin=0;ixpr=0;mxstep=0;maxordn=0;maxords=0;
  447.  
  448. //provisional values as default
  449.  
  450. jacflag=2;ml=-1;mu=-1;
  451.  
  452. %ODEOPTIONS=[itask,tcrit,h0,hmax,hmin,jacflag,mxstep,maxordn,maxords,ixpr,ml,mu];
  453.  
  454. %ODEOPTIONS(7)=10;
  455.  
  456. //ode(x0,t0,t,f);  // ---> Non convergence
  457.  
  458.  
  459. //            Test of %ODEOPTIONS(8:9)
  460.  
  461. //%ODEOPTIONS(8:9)=[maxordn,maxords] ---> maximum order for nonstiff and stiff
  462.  
  463. itask=1;tcrit=0;h0=0;hmax=0;hmin=0;ixpr=0;mxstep=0;maxordn=0;maxords=0;
  464.  
  465. //provisional values as default
  466.  
  467. jacflag=2;ml=-1;mu=-1;
  468.  
  469. %ODEOPTIONS=[itask,tcrit,h0,hmax,hmin,jacflag,mxstep,maxordn,maxords,ixpr,ml,mu];
  470.  
  471. %ODEOPTIONS(8:9)=[0,0]; //---> default values
  472.  
  473. wref=ode(x0,t0,t,f); //just for computing reference
  474.  
  475.  
  476. %ODEOPTIONS(8:9)=[4,3];
  477.  
  478. ww=ode(x0,t0,t,f);norm(wref-ww,1);
  479.  
  480.  
  481. %ODEOPTIONS(8:9)=[12,5];
  482.  
  483. if norm(wref-ode(x0,t0,t,f),1)>Leps then bugmes();quit;end
  484.  
  485.  
  486. //using stiff method
  487.  
  488.  
  489. %ODEOPTIONS(9)=0;
  490.  
  491. wref=ode('st',x0,t0,t,f);
  492.  
  493. %ODEOPTIONS(9)=5;
  494.  
  495. if norm(wref-ode('st',x0,t0,t,f),1) >Leps then bugmes();quit;end //=0
  496.  
  497. %ODEOPTIONS(9)=4;
  498.  
  499. if norm(wref-ode('st',x0,t0,t,f),1)  >Leps then bugmes();quit;end  //small
  500.  
  501.  
  502.  
  503. //using nonstiff method
  504.  
  505.  
  506. %ODEOPTIONS(8)=0;
  507.  
  508. wref=ode('ad',x0,t0,t,f);
  509.  
  510. %ODEOPTIONS(8)=12;
  511.  
  512. if norm(wref-ode('ad',x0,t0,t,f),1) >Leps then bugmes();quit;end  //=0
  513.  
  514. %ODEOPTIONS(8)=5;
  515.  
  516. if norm(wref-ode('ad',x0,t0,t,f),1) >Leps then bugmes();quit;end   //small
  517.  
  518.  
  519. //mixed
  520.  
  521. %ODEOPTIONS(8:9)=[5,12];
  522.  
  523. wref=ode(x0,t0,t,f);
  524.  
  525. %ODEOPTIONS(8:9)=[4,10];
  526.  
  527. if norm(ode(x0,t0,t,f)-wref,1)>Leps then bugmes();quit;end   //small
  528.  
  529.  
  530.  
  531. A=diag([-10,-0.01,-1]);
  532.  
  533. deff('uu=u(t)','uu=sin(t)');
  534.  
  535. B=rand(3,1);
  536.  
  537. deff('y=f(t,x)','y=A*x+B*u(t)')
  538. Warning :redefining function: f                       
  539.  
  540.  
  541. %ODEOPTIONS(1)=2;
  542.  
  543. yy1=ode('stiff',[1;1;1],0,1,f);
  544.  
  545. yy2=ode('stiff',[1;1;1],0,2,f);
  546.  
  547. %ODEOPTIONS(1)=3;
  548.  
  549. yy1=ode('stiff',[1;1;1],0,1,f);
  550.  
  551. yy2=ode('stiff',[1;1;1],0,2,f);
  552.  
  553.  
  554. clear %ODEOPTIONS;
  555.  
  556. rand('seed',0);rand('normal');
  557.  
  558. nx=20;A=rand(nx,nx);A=A-4.5*eye;
  559.  
  560. deff('y=f(t,x)','y=A*x')
  561. Warning :redefining function: f                       
  562.  
  563.  
  564. deff('J=j(t,x)','J=A')
  565. Warning :redefining function: j                       
  566.  
  567.  
  568. //%ODEOPTIONS(1)=1;
  569.  
  570. y2=ode('stiff',ones(nx,1),0,2,f,j);
  571.  
  572. [y1,w,iw]=ode('stiff',ones(nx,1),0,1,f,j);
  573.  
  574. y2p=ode('stiff',y1,1,2,f,j,w,iw);
  575.  
  576. y12=ode('stiff',ones(nx,1),0,[1,2],f,j);
  577.  
  578. norm(y12(:,2)-y2p);
  579.  
  580. yaf=ode('adams',ones(nx,1),0,2,f,j);
  581.  
  582. yaj=ode('adams',ones(nx,1),0,2,f,j);
  583.  
  584. ysf=ode('stiff',ones(nx,1),0,2,f,j);
  585.  
  586. ysj=ode('stiff',ones(nx,1),0,2,f,j);
  587.  
  588.  
  589.  
  590.